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Abstract: A method is detailed for the phase space integration of 7* —> multi-jets ap¬ 
plicable to parton level Monte Carlo’s at any order in perturbation theory. Other non-jet 
objects, massless or massive, can be included in the phase space generation. We corre¬ 
late the bremsstrahlung events in a manner that integrates out all partonic configurations 
leading to a fixed jet configuration, thereby improving convergence. This also allows the 
method to extend infra-red safety to the fully differential multi-jet cross section. 
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1. Introduction 

When calculating differential cross sections involving jets the majority of effort and in¬ 
genuity is invested in evaluating the matrix elements. The phase space integration and 
definitions of proper observables are usually afterthoughts. For Leading Order (LO) par- 
ton level Monte Carlo’s many methods have been proposed to improve the behavior of 
multi-jet event generation. These methods focus on the properties of the LO amplitude [1- 
7]. Here we focus on event generation for parton level event generators beyond LO. At 
higher order partons are clustered, resulting in jets composed of two or more partons. This 
internal configuration of the jet is not purview to the observer within the context of a 
fixed order calculation and needs to be averaged, resulting in a perturbative calculable and 
reliable result. 

We will exploit this necessity by constructing a phase space generator which explicitly 
integrates out all partonic configurations which give rise to the same multi-jet finals state 
for a particular jet algorithm. In ref. [8] the basic concepts of such an approach to phase 
space integration was developed and a proof of existence was given. These techniques were 
refined and modified in ref. [9] for application to the Matrix Element Method (MEM), 
extending the concept of MEM to Next-to-Leading Order (NLO) [10]. Finally, a first look 
at how to implement these correlated phase space generators for the calculation of NLO 
cross sections were highlighted in ref. [11]. In this paper we develop the method fully for 
parton event generators at lepton colliders, ready to be applied to e.g. fully differential 
Next-to-Next-to-Leading Order (NNLO) 3-jet cross sections [12]. 

The standard manner in which phase space event generators work in Monte Carlo’s cal¬ 
culating higher order correction to jet cross sections is to generate the events with different 
parton multiplicities completely uncorrelated. For example, to calculate the NNLO 3-jet 
cross section at a leptonic collider one first generates many 3-parton events and bins the ob¬ 
servable calculated from the 3 momenta, weighted by the value of the sum of the regulated 
3-parton amplitudes. Once completed, another run is initiated generating many 4-parton 
events for which the observable is calculated and binned with the weight of the regulated 
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4-parton amplitude. Finally, the 5-parton events are generated and the observable binned 
using the regulated 5-parton amplitude. In this the integration of the observable over the 
jet phase space plays a crucial role to achieve infra-red safety. 

Yet, from the viewpoint of the observed jet final state, these different parton multiplic¬ 
ity events have a high degree of correlation. By applying a jet algorithm, the hadronic final 
state is simplified into a final state of jet objects, each jet having a momentum. Perturba¬ 
tive QCD should be able to calculate these jet cross sections, i.e. the correlations between 
the jet momenta, order by order in perturbation theory provided the jets are opaque. That 
is, the jets are averaged over all hadronic configurations leading to that particular jet con¬ 
figuration. As a consequence the fully differential jet cross section should be calculable 
given an appropriate infra-red safe jet algorithm. For the NNLO 3-jet example this will 
significantly alter the manner in which the 3-parton, 4-parton and 5-parton final states 
are generated. First of all, the generator calculates the fully differential 3-jet cross section 
d < 73 / dpidp- 2 <ipz where {pi} are the jet axis momenta. Observables can be calculated by 
integrating the jet observable over the jet phase space. Note that nowhere it is required 
to average over the jet phase space to obtain an infra-red safe answer, nor are there any 
constraints on the type of jet observable one can look at. For the event generator this 
means that the starting point is a given fixed 3-jet configuration (i.e. the jet momenta 
{pi} are given). This requires only a single evaluation of the 3-parton amplitudes. From 
this starting point we generate 4-parton events which reconstruct back to the initiating 
3-jet momenta using the jet algorithm. Because the jets are fixed the 4-parton amplitude 
weights can simply be added to the 3-parton amplitude weight. Due to the recursive nature 
of sequential clustering in jet algorithm we can repeat the above argument for the 5-parton 
contribution: given a single 4-parton event we can generate many 5-parton events such 
that applying a single step in the jet algorithm leads back to the initiating 4-parton event. 
It is clear that this will lead to a highly correlated event generator and the event weight of 
the fully differential 3-jet cross section is a series expansion in the strong coupling constant 
at the scale of the jet resolution. 

Current jet algorithms do not extend the concept of infra-red safety to the fully dif¬ 
ferential jet cross sections. This is solely due to the clustering phase in the algorithm, 
i.e. how to combine two momenta to form a new one. The new momentum formed from 
the momenta of the respective clustered particles is simply the sum of these momenta. 
For some algorithms this new momentum is redefined to make it a massless momentum 
(see ref. [13] for an overview of these schemes). As a result applying the jet algorithm 
on e.g. a 4-parton final state to form a 3-jet final state will never overlap with the 3-jet 
final state generated from the 3-parton final state due to jet mass and/or non-momentum 
conservation. 

As a consequence we have to adjust the cluster phase of the jet algorithms in order 
to obtain infra-red safe fully differential jet cross sections. Note that one can still apply 
any ordinary jet algorithm using the above described correlated phase space generator. 
However the now theoretical jet algorithm used internally by the generator to reorganize 
the correlated phase space generation does not match the applied jet algorithm exactly. 
The price to pay is that one must define observables and average this jet observable by 


- 2 - 



integrating over the full jet phase space so that an infra-red safe prediction can be made. 
Still, the generation of the correlated multi-parton final states could greatly benefit the 
convergence of the predictions of the observable made by the Monte Carlo. 

In section 2 we will construct the required forward brancher which converts an n- 
particle phase space into a (n + l)-particle phase space, invertible using the cluster al¬ 
gorithm. In section 3 this Forward Branching Phase Space (FBPS) generator is used to 
construct the correlated event generators needed for calculating fully differential multi-jet 
cross sections at any order in perturbation theory. We conclude in section 3 by summarizing 
the results and outline the next steps in further developing the FBPS method. 

2. Constructing the forward branching phase space generator 

The first step in the construction of the FBPS is to define a proper clustering algorithm. 
As already explained in the introduction, we need a clustering algorithm which maps a 
massless n-particle phase space onto a massless (n — l)-particles phase. 1 

Current jet algorithms simply combine two momenta by adding the 4-vectors. Take as 
an example the decay of a heavy particle with momentum Q into four massless particles, 
clustering momentum p 3 and momentum p± 

Q =Pi+P2+P3+P4=Pi+P2+P34 ■ (2.1) 

This obviously lead to a 3-particle phase space where one of the momenta is massive. Some 
clustering algorithms rescale the energy or momentum such that the clustered particle is 
massless (see ref. [13]). However such a scheme would change Q through momentum 
conservation. 

The way to modify the clustering is to use a recoil momentum. By rescaling at least 
one of the other non-clustered particles one can make the clustered particle massless (this 
is in essence a 3 —> 2 clustering as is commonly used to construct the subtraction terms 
needed to regulate various amplitudes in higher order calculations [14]). Explicitly, eq. 2.1 
becomes 

Q = P1+P2+P3+P4 = P1+P2+P34 = {1-/3)(P1+P2) + (P34:+P(P1+P2)) = P1+P2+P3 , (2.2) 

where we have used the sum of all other particles as a recoil momentum so to avoid further 
complicating the algorithm on how to select a recoil particle. The new momenta are given 
by pi = (1 — f3)p\, P 2 = (1 — (3)p2 and p 3 = p^ + @(pi + $ 2 )- The branching scale variable 
/3 can now be chosen such that the clustered particle is massless 

P\ = (P34 + PP\2) 2 = P 2 S\2 + 2/3 P 34 ' Pl2 + p\i = 0 

(P34 ■ P12) ~ J(P34 ■ P12) 2 ~ Sl2Pu 
—^ R — v 


1 Note that it is straightforward to add other, massive or massless, momenta as long as they to not 

participate in the jet clustering. For the remainder of the paper we ignore this in order to simplify the 
notation and discussion. 


(2.3) 
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were pij = pi + Pj and Sij = pf , J . With this simple augmentation of the cluster algorithm 
we have what we need to construct the cluster invertible FBPS generators. 

It is convenient to introduce here the inverse of the clustering kinematics, i.e. the 
branching kinematics. Introducing the branching scaling variable a we have 


Q = P1+P2+P3 = (l + a)(pi+p 2 ) + (p3~a(pi+p2 )) = P1+P2+P34 = P1+P2+P3+P4 , (2.4) 
where a is given by the quadratic equation 

«34 = P34 = (p 3 - a(pi + p 2 )) 2 = a 2 s\2 - 2ap 3 ■ {p\ + p 2 ) • (2.5) 

The branching and merging scale variables are related by 


1-/3 = 


1 


1 OL 


P = 


a 


1 + a ’ 


and give the relations to be used later 


13{pi +p 2 ) = f3p\2 = /3( 1 + a)pi 2 = apu = a{pi + p 2 )] d/3 = 


da 


(1 + a ) 2 


( 2 . 6 ) 


(2.7) 


Having augmented the clustering algorithm we can now start making the FBPS gener¬ 
ator needed to generate the massless (n-fl)-parton phase space from the massless n-parton 
phase space subject to the clustering constraint of eq. 2.2. 

The n-particle fully differential cross section is given by 


da 


( (2vr) 4 \ 


x d$ n (Q;pi,... ,p n ) 


dpi---dp n yiy/Q 1 ) 

The flat phase space is given by 

n d A n- 

d<& n {Q\PU---,Pn) = Ilj^- 3 S + ( P 1)S\Q-Pl - Pn) 


( 2 . 8 ) 


(2.9) 


We will derive the FBPS generating a massless 4-particle phase space from a massless 
3-particle phase space. Afterward we will generalize this to n-particle phase spaces. The 
first step in constructing the FBPS generator is simply generating the massless 4-particle 
phase space from a 3-particle phase space with one massive particle 


d ( I>4(<3;pi P2P3IM) = (27r) 3 ds 3 4 4>3(Q;PlP2P34) x ^ 2 {P34',PlP2) 


( 2 . 10 ) 


where p| 4 = S34. 

The next step is to implement the rescaling of momentum P34 using the recoil momen¬ 
tum pi 2 - Specifically, starting from the massive phase space 

d4> 3 (Q;pip 2 P34) = ^^^^^pd + (pf)(5 + (^)5 + (p| 4 -S34)d 4 (Q-pi-p2-P34) , (2.11) 
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we have to derive the Jacobean J generated by the change of the integration momenta 


d^ 3 {Q-pip 2 pu) = J x d$ 3 (Q-,pip 2 p 3 ) (2.12) 

= JX ^3 ( 2^3 ( 2 ^3 S+ (Pi) 6+ (P 2 2) S+ (pI) 64 (Q -P1-P2- P3) ■ 

To calculate the Jacobian J we need to rewrite the integrals over the momenta pi,p 2 
and P 34 into integrals over momenta pi, p 2 and p 3 . To do this we mathematically express 
the clustering algorithm as a decomposition of unity 

1 = S 12 (/3+-/3_) J dpdpidp 2 dp35((l- (})pi-pi)6((l- P)p 2 -p2)S(p3) 

X d(p 3 - (p 34 + PP12)) , (2.13) 

where /3± are given by solving equation {p 3 4 + fdpu) 2 = 0 : 

si2((3+ ~ P-) = 2y/(p34 • P 12) 2 - S 34 S 12 • (2-14) 


The above identity is expressed in terms of the cluster scale variable (3 and cluster kine¬ 
matics. We rewrite this equation in terms of the branching scale variable a and branching 
kinematics, giving us 

1 = 2 J da dpi dp 2 dp 3 {l + a)'\/(P 34 -pi 2 ) 2 - S 34 M 2 

x 6{pi - (1 + a)pi)5{p2 - (1 + oc)p 2 )6{pl)S{pu - (p 3 - ap i2 )) ■ (2.15) 

Multiplying this decomposition of one to eq. 2.11 and integrating the momenta p \, p 2 and 
P 34 over the appropriate 5-functions one obtains 


d$3(Q',PiP2P34) = J da{ 1 + a) 3 i/(p 3 4 ■ P 12) 2 - S 34 S 12 5(a 2 si2 - 2 ap 3 ■ p i2 - s 34 ) 

dpi dp 2 dp 3 


(27t) 3 (27t) 3 (27t) 3 


d{Pi)d{p 2 )d{p 3 )S{Q -P 1 -P 2 - P3) 


(2.16) 


with pi = (1 + a)pi, P 2 = (1 + a)p 2 and j5 34 = p 3 — api 2 . 

By integrating over the branching variable a and selecting the physical solution of the 
quadratic equation we obtain the final expression 


d^ 3 {Q;pip2Pu) 


with 


n , 'iS (P34-P12) -S 34 

(1+<L) V(ft»w+»» 


x d^ 3 {Q\pip 2 p 3 ) , 


(2.17) 


P34 • P12 — V (^34 • P12) 2 + S34S12 

a- =-; pi = (l+a_)pi; p 2 = (l+a_)p 2 ; P34 =P3~oi-Pi2 ■ 

S 12 

(2.18) 

It is now straightforward to generalize this result. We can generate the massless 71- 
particle phase space from a massless (n — l)-particle phase space using the FBPS generator. 
Specifically, by branching massless particle j we obtain 

d<S> n (Q-,{p}) = d<S> n -i{Q;{p}) x , (2.19) 
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where 2 


d$ Sps((P}l(P}) 


Pi 

Pjn 

Qj 


Oi 7 


(2vr) a (l + aj ) 2n - 3 x 


( Pjn'Qj ) Sj n Qj , A A 

" ds jn Cl < I ) 2 \Pjn ; PjPr 


(Pj ' Qj) “b Sj n Q 


(1 + a j)pi (i^j) 

Pj - ajQj = (1 + otj)pj - ajQ 

Ek = Q~Pj 


i+j 


(pj ■ Qj) \J(pj • Qj) + SjnQ 


Q 2 j 


( 2 . 20 ) 


We introduced the notation {p} = Pi,P 2 , ■ ■ ■ for a list of momenta in order to simplify the 
notation. The number of momenta in the list is clear from the context. 

By iterating the above FBPS, we can generate a n-particle phase space from a m- 
particle phase space (with at least 1 massless particle) through intermediary fc-particle 
phase spaces. Furthermore, we can invert the generated n-particle phase space by applying 
the appropriate clusterings as defined in eq. 2.2. 




Figure 1: Comparison of the flat phase space distribution and the 3 —> 4 FBPS generated dis¬ 
tribution for the observable ma Xi(Ei)/y/Q^ using 10 8 events. The figure on the right is the ratio 
of the two distributions. Both the 3-particle and 4-particle flat phase space were generated using 
RAMBO and \/Cp = 240. The results in the left graph were not rescaled for the bin-width. 

Finally we want to validate the correctness of eq. 2.19. This can be done using a Monte 
Carlo integration of the FBPS generated results versus the results obtained using the flat 
phase space generator RAMBO [15]. 

The first check is to see whether the FBPS generated phase space has the correct phase 


2 If we start from the special configuration d&(Q,pip 2 ) where both momenta are massless, we get the 

solution a = — S23/S12 = —S 23 /Q 2 when branching momentum p 2 - 
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Figure 2: Comparison of the flat phase space distribution generated by and the 3 —► 4 —>• 5 
generated FBPS distribution for the observable maXj(Fj )/yfcft using 10 s events. The figure on the 
right is the ratio of the two distributions. Both the 3-particle and 5-particle flat phase space were 
generated using RAMBO and y/Cp = 240. The results in the left graph were not rescaled for the 
bin-width. 


space volume. The volume of phase space is given by 


f - /d<M<2;{p (m) }) X 


n ({p (ro+l) }io> (m+ '=- 1) }) 


L k=\ 


/yrxn-l 1 ( Q 2 ) n 2 

\2/ (27r) 3ri (n — l)!(n — 2)! 


( 2 . 21 ) 


We numerically verified the thus generated n-particle phase space, for all 3 < m < n < 10, 
gives back the correct phase space volume. For each event generated in the Monte Carlo 
all jk were chosen randomly. The m-particle phase space in the comparison was generated 
using the flat phase space generator RAMBO. 

To further validate the FBPS generator we made the differential cross section for the 
observable maxj(Fij)/ y/Q^ where Ej is the energy of momentum pi for both the 3 —» 4 FBPS 
generator in figure 1 and the 3 —> 5 FBPS generator in figure 2. The results were compared 
against the predictions from RAMBO for the 4- and 5-particle phase space respectively. 
As can be seen in both cases the agreement is as expected, given the 10 8 generated events 
used in the Monte Carlo’s. 


3. Constructing differential jet cross sections 

Now that we have the basic building block in the form of the FBPS generator as detailed 
in eq. 2.19, we can construct the event generator to calculate the fully differential jet cross 
sections. The event generator will depend on the resolution function dij = d(pi,pj) of 
the jet algorithm. The pair of momenta with the smallest djj will be clustered by the 
jet algorithm and it will keep clustering pairs of momenta until it reaches the condition 
min ij(dij) > d CVLt (exclusive jet cross section) or until a certain multiplicity of momenta is 
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reached (i.e.. keep clustering until we reach m momenta, giving us an inclusive m-jet cross 
section). 

We will start with eq. 2.19 and extend it to give us the generator to calculate the 
fully differential jet cross sections. The jet resolution function will do two things for 
us. It will give which momentum to branch and set an upper integration boundary on 
the variable Sj n and decay products of pj n . We can implement the branches by using the 
resolution function to partition the (n — l)-particle phase space into wedges. Each wedge 
is associated with a momentum to branch in order to generate the n-particle phase space. 
To do this we define the following partition of one 

n —1 

1 = ^ 6(d jn = d„un) , (3.1) 

3 = 1 

where the 0-function equals one if its argument is true and zero otherwise. The argument 
of the 0-function d m \n = min ij{dij) = min ij(d(pi,pj)) is true as long as dj n has the smallest 
resolution parameter (i.e. it will be the pair which will be clustered by the jet algorithm). 
This partition divides phase space in the desired wedges. In wedge j, particle pj branches 
to give particles {pj,p n } in such a manner that when applying the jet algorithm it will pick 
this pair to cluster back to pj as dj n has the smallest resolution parameter. Note that this 
partitioning include n-jet final states if we demand a clustering cutoff when d m i n > d c u t- 
This distinction can be implemented as 

n— 1 

1 = 'y ( @(djn = d m in) ^0(dmin ^ d c ut) T 0(dmin ^ d C ut)^ • (3-2) 

3 = 1 

The (n — l)-particle jet exclusive FBPS generator is now given by multiplying eq. 2.19 with 
the above partition of one 

d<S>n(Q]{p}) = d$„_i(Q;{p}) x d$fops({£}|M) (3.3) 


d$frp S ({p}\{p}) = Y = 0min)0(dmin < dcut)d $£ ps ({p}|{p}) ■ (3-4) 

3 =1 

This FBPS generator will produce all bremsstrahlung radiation momenta for a fixed (n— 1)- 
jet configuration with jet momenta {pi} provided the modified clustering of eq. 2.2 is used. 
By integrating over the jet phase space and removing the d cut -constraint, the full n-particle 
phase space is obtained. We can iteratively apply this FBPS generator to obtain the 
multiple bremsstrahlung phase spaces. 

To define the fully differential jet cross section at NNLO using the improved clustering 
is now straightforward. Given the n-jet configuration with massless jet momenta {p} we 
get 


NNLO 

a u impr 

d 3 pi ■ ■ ■ d 3 p n 


-4({p})+ / rf®Sg({p}IW) x 


S({p})+ / <if£pi({#}|{ji})C({ji}) 


(3.5) 
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where 


A({ P }) 

cm) 


A{pi ■ • • Pn ) = M {0 \pi ■ • • Pn ) +X (1) (pi • • -Pn) + X (2) (pi ' ' -£>n) 
S(Pl--'Pn) = X ( 0 ) (Pl-"Pn)+X ( 1 ) (Pl"-Pn) 


C(Pl---Pn) = 




l'-'Pn 


(3.6) 


are the appropriate set of matrix elements at each parton multiplicty. Note that d <3?^' is 
a 3-dimensional integral, so that the above integral is at most 6-dimensional. 

We can still use the FBPS generator and apply other jet algorithms to its produced 
events. There should still be an advantage over the usual uncorrelated generation of the 
different phase spaces. Given the jet algorithm mapping A({p}|{]5}) the differential cross 
section is given by 


d 3 V£ NLO 

d 3 pi ■ ■ ■ d 3 p n 


^(M)A(WIM) 


+ /d<i>f p c i({p}|M)A(M|{p}) X 


B({p}) + 


d^imm)^p}m)cm) . 

(3.7) 


However this fully differential cross section is ill-defined and not infra-red safe. Suppose we 
choose all the jet momenta massless, then the bremsstrahlung events will never contribute 
because the bremsstrahlung events generate either jets with masses or momentum Q is 
changed. That is, in this case the virtual and bremsstrahlung events do not merge. To 
obtain infrared safety in this case we must define an appropriate observable such that the 
bremsstrahlung and virtual contributions are sufficiently merged, that is 


dal NL ° 

dO 


r rl 3n rr NNLO 

J d Pi ‘ ’ d 


(3.8) 


In this sense the improved jet algorithm is superior to the standard jet algorithms as it 
extends infra-red safety to the fully exclusive jet cross section. 

We can generalize eq. 3.5 and 3.7 using a generating functional T 


d 3r 


o\ 


impr 


d 3 pi ■ ■ ■ d 3 p r 


— wd({p}) x Ti mpr ({p}); 


d 3r 


O’A 


d 3 pi ■ ■ ■ d 3 p r 


= -A({P}) x r A(W) 


(3.9) 


to generate the n-th order fully differential cross section calculator by expanding n times 
the T-functional 


iWIM) = i + / diCWIIM) x 
ra(W) = 1 +f rf*Sg({#}IW) X A({p}|{p}) x (^|y) • (3.10) 


These generators can readily be implemented into Monte Carlo programs. 
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4. Conclusions 


We constructed a correlated bremsstrahlung parton level event generator which should lead 
to faster convergence of Monte Carlo phase space integrations for higher order jet cross 
sections at lepton colliders. The bremsstrahlung events are generated highly correlated, 
enabling better cancellations between real and virtual contributions. One can use any jet 
algorithm on the generated events, however using the augmented jet algorithm enables us 
to define the fully exclusive multi-jet differential cross section. That is, for any given jet 
configuration the generator will integrate out all radiation inside the jets, thereby mak¬ 
ing them opaque, and combine the matrix element weights of the different multiplicities 
unencumbered by any constraint. In other words we can define a probability density (or 
K -factor), calculable order by order in perturbation theory, for every exclusive jet config¬ 
uration. 

The next step is to extend this method to hadron colliders and apply the above method 
to NLO and NNLO multi-jet generators at lepton and hadron colliders. 
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